We are IntechOpen, the world's leading publisher of Open Access books Built by scientists, for scientists

Open access books available 5,300

130,000 155M

International authors and editors

Downloads

Our authors are among the

most cited scientists 154 TOP 1%

Selection of our books indexed in the Book Citation Index in Web of Science™ Core Collection (BKCI)

# Interested in publishing with us? Contact book.department@intechopen.com

Numbers displayed above are based on latest data collected. For more information visit www.intechopen.com

# **Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals**

Roman Antos <sup>1</sup> and Martin Veis 1,2 1 *Institute of Physics, Faculty of Mathematics and Physics, Charles University in Prague* 2 *Institute of Biophysics and Informatics, First Faculty of Medicine, Charles University in Prague Czech Republic*

# **1. Introduction**

Photonic crystals are modern artificially designed periodical systems capable to affect the motion of photons in a similar way that the periodicity of the atomic potential in a semiconductor crystal affects the motion of electrons. The physical properties of light in a photonic crystal resemble those of electrons in atomic crystals, leading to forbidden propagation of electromagnetic modes at certain frequencies, as demonstrated by Yablonovitch (1987), Ho et al. (1990), Joannopoulos et al. (1995; 1997), etc. The existence of the optical band gap (which is the part of the spectrum for which the wave propagation is not possible) makes photonic crystals broadly interesting from many viewpoints of fundamental research and applications. Recent studies are motivated by promising applications such as purely optical integrated circuits [e.g., by Lin et al. (1998), Noda (2006), or Hugonin et al. (2007)], artificial metamaterials with high tunability [Datta et al. (1993); Genereux et al. (2001); Krokhin et al. (2002); Reyes et al. (2005)], high-sensitivity photonic biosensors [Block et al. (2008); Skivesen et al. (2007)], or devices based on phenomena not accessible in conventional media [Benisty (2009); Kosaka et al. (1998); Krokhin & Reyes (2004)]. It is also necessary for accounting for structural colors of wings of butterflies or beetles, feathers of birds, or iridescent plants [Kinoshita & Yoshioka (2005); Vukusic & Sambles (2003)].

Since the optical properties of photonic crystals strongly depend on their geometrical structure, used materials, etc., their proper design is crucial for the correct device functionality. Thorough theoretical analysis therefore takes place in the development, using various numerical methods for calculation including methods of finite difference in the time or frequency domain or finite element methods [Joannopoulos et al. (1995)]. One of the most common calculation techniques applied to photonic crystals is the plane wave expansion method, which is a frequency-domain approach based on the expansion of the fields and material parameters into the Fourier (or reciprocal) space. The components of this expansion represent the definite-frequency states. After some necessary truncation of the complete basis (plane waves with a finite cutoff), the partial differential equations are then solved as a linear-algebraic problem. However, the convergence rate of this method strongly depends on the implementation of Maxwell's equations in the truncated plane-wave basis [Meade et al. (1993); Sozuer et al. (1992)]. In the case of periodic discontinuities (typical for photonic crystals) the convergence is rather poor so that the computer calculations might become extremely time- and memory-consuming.

Because the underlying physical phenomenon in the optical behavior of photonic crystals is based on diffraction (therefore the lattice constant of a periodic structure has to be in the same length-scale as half the wavelength of the electromagnetic wave), several conclusions can be advantageously adopted from the classical coupled-wave theory, one of the most effective methods for modeling diffraction of electromagnetic waves by periodic gratings, which was developed during several past decades. This can provide a high enhancement to the plane wave expansion method, resulting in the reduction of the computation resources.

In the mid 1990s Li (1996) showed that Laurent's "direct" rule, which had always been adopted in conventional formulations to factorize the truncated Fourier series that corresponds to products of two periodic functions, presents bad convergence when the two functions of the product are simultaneously discontinuous. He suggested three Fourier factorization rules (briefly summarized in Section 2.3) and applied them to one-dimensional (1D) diffraction gratings. This major breakthrough in the grating theory (called "fast Fourier factorization") was soon applied by many authors to various grating structures with arbitrary periodic reliefs, anisotropic [Li (1998)] and slanted [Chernov et al. (2001)] periodic systems, their various combinations [Li (2003); Watanabe (2002); Watanabe et al. (2002)], and other systems [Bonod et al. (2005a;b); Boyer et al. (2004)].

Later Li (1997) applied the factorization rules to two-dimensional (2D) periodic structures treated by "zigzag" Fourier expansion, which yielded an improvement for rectangular dots or holes. However, Popov & Neviere (2000) have pointed out that the staircase approximation (of the coupled wave theory using the slicing of relief profiles) in combination with the traditional formulation of differential equation within one slice violates Li's factorization rules. This was a major complication for the analysis of the periodic systems made of rounded elements. Therefore, they applied a coordinate transform to treat individually the normal and tangential components of the electric field on 1D sinusoidal-relief gratings, which enabled the application of the correct rule for each field component and thus improved the convergence.

Later David et al. (2006) utilized the normal–tangential field separation to 2D photonic crystals composed of circular or elliptical holes. Similarly, Schuster et al. (2007) applied this method to 2D gratings, and also suggested more general distributions of polarization bases [Gotz et al. (2008)]. These approaches, always dealing with linear polarizations, enabled a significant improvement of the convergence properties, but ignored the fact that the transformation matrix between the Cartesian and the normal–tangential component bases of polarization became discontinuous at the center and along the boundaries of the periodic cell, which slowed down the resulting convergence. To overcome these discontinuities, a distribution of more complex (i.e., generally elliptic) polarization bases was recently suggested to improve optical simulations of 2D gratings and photonic crystals [Antos (2009); Antos & Veis (2010)].

Our chapter will describe in detail the application of the Fourier factorization rules to the plane wave expansion method for numerical analysis of general photonic crystals. Section 2 will introduce the principle of the plane wave expansion together with the notation of matrices and factorization theorems. Section 3 will refer to 1D photonic structures made as periodic stratified media. The consistency of the correct factorization rules with classical theory of Yeh et al. (1977) and Yariv & Yeh (1977) will be shown, pointing to the correct boundary conditions of the tangential components of the electric and magnetic field on multilayer interfaces. Section 4 will repeat our previously described methodology for 2D photonic crystals made of circular elements, and Section 5 will generalize it to elements of other shapes. Sections 6

and 7 will propose how to factorize anisotropic and three-dimensional (3D) photonic crystals, respectively.

### **2. Plane wave expansion**

#### **2.1 General remarks**

The modes of photonic crystals are in principle the eigensolutions of the wave equation with an inhomogeneous, periodic relative permittivity *ε*(r). One possible version of the wave equation is the equation for the unknown electric field with an unknown frequency,

$$\stackrel{\smile}{\frac{1}{c^2}}\nabla \times (\nabla \times \mathbf{E}) = \frac{\omega^2}{c^2} \mathbf{E'} \quad \stackrel{\smile}{\quad} \stackrel{\smile}{\quad} \stackrel{\smile}{\quad} \stackrel{\smile}{\quad} \quad \rightarrow \tag{1}$$

where *ω*2/*c* 2 is its eigenvalue (*ω* is the frequency and *c* the light velocity in vacuum) and

$$E(r) = e^{-ik \cdot r} \sum\_{m,n,l} e\_{mnl} e^{-ik\_0(mpx+nqy+lsz)} = \sum\_{m,n,l} e\_{mnl} e^{-ik\_0(p\_mx+q\_ny+s\_lz)} \tag{2}$$

is its eigenfunction, which has the form of a pseudoperiodic Floquet–Bloch function. Here *p* = 2*π*/Λ*x*, *q* = 2*π*/Λ*y*, and *s* = 2*π*/Λ*<sup>z</sup>* are the normalized reciprocal lattice vectors. For simplicity we assume the periods Λ*<sup>j</sup>* along the Cartesian axes throughout this chapter. For brevity we have also defined k = *k*0[*p*0, *q*0,*s*0], *p<sup>m</sup>* = *p*<sup>0</sup> + *mp*, *q<sup>n</sup>* = *q*<sup>0</sup> + *nq*, and *s<sup>l</sup>* = *s*<sup>0</sup> + *ls*. (Analogously we could write the wave equation for an unknown magnetic field H or any other field from Maxwell's equations.)

Owing to the periodicity of the problem, the plane wave expansion method is the reference method for the mode calculation. It is based on the Fourier expansion of the field such as in Equation 2 and on the Fourier expansion of a material function, either the permittivity or the impermittivity *η*(r)= 1/*ε*(r),

$$\varepsilon(\mathbf{r}) = \sum\_{m,n,l} \varepsilon\_{mnl} e^{-ik\_0(mp\mathbf{x} + nqy + l\mathbf{s}z)} \tag{3}$$

$$\eta(\mathbf{r}) = \sum\_{m,n,l} \eta\_{mnl} e^{-ik\_0(mp\mathbf{x} + nqy + l\mathbf{s}z)} \tag{4}$$

The rules for choosing the most appropriate material parameter and the most appropriate field for the Fourier expansion are governed by various methods of Fourier factorization. In the past, the E method (*η* and the electric displacement D were expanded), H method (*η* and H were expanded), and Ho method (*ε* and E were expanded) were the typical choices.

#### **2.2 Matrix notation**

Now we carry out the transformation of the partial differential equations into matrix equations in order to solve the eigenproblem by linear-algebraic methods. For simplicity we limit ourselves to 1D and 2D photonic crystals, and always choose the direction of propagation in the *xy* plane, so that *∂<sup>z</sup>* = 0. With these restrictions we write

$$\varepsilon(x,y) = \sum\_{m,n=-\infty}^{+\infty} \varepsilon\_{mn} e^{-i(mpx+nqy)}\,,\tag{5}$$

$$f(x,y) = \sum\_{m,n=-\infty}^{+\infty} f\_{mn} e^{-i(p\_m x + q\_n y)}.\tag{6}$$

where *f* is a component of the electric field. We will now derive the matrix expressions of the fundamental relations

$$h(\mathbf{x}, y) = \varepsilon(\mathbf{x}, y) f(\mathbf{x}, y),\tag{7}$$

$$g\_{\chi}(\mathbf{x}, y) = \partial\_{\chi} f(\mathbf{x}, y), \tag{8}$$

$$\mathbf{g}\_{\mathcal{Y}}(\mathbf{x}, \mathbf{y}) = \partial\_{\mathcal{Y}} f(\mathbf{x}, \mathbf{y})\_{\prime} \tag{9}$$

i.e., the relations of multiplication by a function and applying partial derivatives. Assuming the expansions of the new functions *h* = ∑*m*,*<sup>n</sup> hmn e* −*i*(*pmx*+*qny*) , *g<sup>x</sup>* = ∑*m*,*<sup>n</sup> gx*,*mn e* −*i*(*pmx*+*qny*) , and *g<sup>y</sup>* = ∑*m*,*<sup>n</sup> gy*,*mn e* −*i*(*pmx*+*qny*) , we rewrite Equations 7–9 using the convolution rule, and applying the partial derivatives as follows:

$$h\_{mn} = \sum\_{k,l=-\infty}^{+\infty} \varepsilon\_{m-k,n-l} f\_{kl} \tag{10}$$

$$g\_{\chi,mn} = -ip\_mf\_{mn} \tag{11}$$

$$g\_{y,mn} = -iq\_nf\_{mn}.\tag{12}$$

Assuming furthermore a finite number of the retained Fourier coefficients, i.e., using the summation ∑ +*M <sup>m</sup>*=−*<sup>M</sup>* <sup>∑</sup> +*N n*=−*N* , we can renumber all the indices to replace the couple of two sets *m* ∈{−*M*, −*M* + 1, . . . , *M*} and *n* ∈{−*N*, −*N* + 1, . . . , *N*} by a single set of indices *α* ∈{1, 2, . . . , *α*max}, with *α*max =(2*M* + 1)(2*N* + 1), related

$$a(m, n) = m + M + 1 + (n + N)(2M + 1),\tag{13}$$

$$n(a) = (a-1)\mathbf{div}(2M+1) - N,\tag{14}$$

$$m(a) = (a-1)\mathbf{mod}(2M+1) - M,\tag{15}$$

where "**div**" denotes the operation of integer division and "**mod**" the remainder (the modulo operation). Then we can rewrite Equations 10–12 into the matrix relations

$$[h] = \|\varepsilon\| [f].\tag{16}$$

$$\begin{bmatrix} \mathbf{g}\_{\mathbf{x}} \end{bmatrix} = -i\mathbf{p} \begin{bmatrix} f \end{bmatrix}\_{\smile \smile} \tag{17}$$

$$\lfloor g\_{\mathbf{y}} \rfloor = \underbrace{-i\mathbf{q}[f]}\_{\frown} \tag{18}$$

where [ *f* ], [*h*], [*gx*], and [*gy*] are column vectors whose *α*th elements are the Fourier [*m*, *n*] elements of the functions *f* , *h*, *gx*, and *gy*, indexed by *α*(*m*, *n*) defined in Equation 13, and where [[*ε*]], **p**, and **q** are matrices whose elements are defined

$$\mathbb{E}[\varepsilon]\_{a\beta} = \varepsilon\_{m(\alpha) - m(\beta), n(\alpha) - n(\beta)}\tag{19}$$

$$p\_{a\beta} = p\_{m(a)} \delta\_{a\beta} \tag{20}$$

$$q\_{\alpha\beta} = q\_{n(\alpha)} \delta\_{\alpha\beta\prime} \tag{21}$$

where the indices on the right hand parts are defined by Equations 14 and 15 and where *δαβ* denotes the Kronecker delta. As a summary we can say that the multiplication by a function is in the reciprocal space represented by the matrix [[*ε*]] (in the sense of the limit *<sup>α</sup>*max → <sup>∞</sup>) and that the partial derivatives are represented by the diagonal matrices −*i***p** and −*i***q**.

For 1D periodicity we choose the inhomogeneity along the *y* axis. In this case the *α* index is not necessary because the *n* index is sufficient,

$$\varepsilon(y) = \sum\_{n=-\infty}^{+\infty} \varepsilon\_n e^{-inqy} \,\_{\prime} \tag{22}$$

$$\begin{array}{ccccc} \begin{array}{c} \Box \\ \Box \end{array} & f(y) = \sum\_{n = -\infty}^{+\infty} f\_n e^{-iq\_ny} & \text{(23)}\\ \Box \end{array} \tag{23}$$

$$\begin{array}{c} \begin{array}{c} \Box \\ \Box \end{array} \end{array} \quad \begin{array}{c} \begin{array}{c} f\_n e^{-iq\_n y} \\ \Box \end{array} \end{array} \quad \begin{array}{c} \begin{array}{c} \Box \\ \Box \end{array} \end{array} \quad \begin{array}{c} \begin{array}{c} \begin{array}{c} \Box \\ \Box \end{array} \end{array} \end{array} \begin{array}{c} \begin{array}{c} \Box \\ \Box \end{array} \end{array} \begin{array}{c} \begin{array}{c} \Box \\ \Box \end{array} \end{array} \begin{array}{c} \begin{array}{c} \Box \\ \Box \end{array} \end{array} \end{array} \tag{24}$$

Here [[*ε*]] is a Toeplitz matrix (a matrix with constant diagonals).

#### **2.3 Simplified theorems of Fourier factorization**

Although the theorems were derived by Li (1996) for 1D periodic functions, we here summarize them in the matrix formalism independent of the number of dimensions. Let *f* , *h*, and *ε* be piecewise-continuous functions with the same periodicity related

$$h = \varepsilon f,\tag{26}$$

and let [ *f* ], [*h*], and [[*ε*]] denote their matrices as defined in Section 2.2.

**Theorem 1.** If *ε* and *f* have no concurrent discontinuities, then the Laurent rule applied to Equation 26 converges uniformly on the whole period and hence

$$[h] = \begin{bmatrix} \varepsilon \end{bmatrix} [f] \tag{27}$$

[ *f* ] (30)

can be applied with fast convergence.

**Theorem 2.** If *ε* and *f* have one or more concurrent discontinuities but *h* is continuous, then Equation 26 can be transformed into the case of Theorem 1,

$$\begin{array}{cccc} & & f = \frac{1}{\varepsilon}h, & & & \text{(28)}\\ \text{and hence} & & & & & \text{(28)}\\ & & & & & \text{(21)}\\ \text{According, we can state} & & & & & & \text{(21)}\\ \text{According, we can state} & & & & & & \text{(21)}\\ & & & & & & \text{(21)} \end{array} \tag{29}$$

referred to as the inverse rule. We say that the functions *ε* and *f* have complementary discontinuities.

*ε*

[*h*]=

**Theorem 3.** If none of the requirements of the first two theorems are satisfied, then none of the rules can be applied correctly because Equations 27 and 30 are no longer valid at the points of discontinuities, which considerably slows down the convergence. Therefore, we should carefully analyze the continuity of the functions and transform all the partial differential formulae to the first two cases.

# **3. One-dimensional photonic crystals**

#### **3.1 Geometry of the problem**

Fig. 1 shows the geometrical configuration of a 1D photonic crystal made as periodic alternation of two different layers with relative permittivities *ε<sup>I</sup>* and *εII* , thicknesses *d*<sup>1</sup> and *d*2, with periodicity Λ along the *y* axis. The relative thickness of the first layer (with respect to the period) is denoted *w*; the relative thickness of the second layer is then 1 − *w*. The coordinate system is chosen to get the uniform problem along the *z* axis. This means that the plane of incidence (plane defined by the vector of periodicity and the wave vector of incidence k = *k*0[*p*0, *q*0,0], here only hypothetical since the photonic crystal is infinite) coincides with the *xy* plane.

Fig. 1. Geometry of a 1D photonic crystal made as periodic alternation of two layers

Then we distinguish between two polarizations of electromagnetic fields which can be treated independently. The transverse electric (TE) polarization has E perpendicular to the plane of incidence (*Ez*, *Hx*, and *H<sup>y</sup>* are nonzero). The transverse magnetic (TM) polarization has H perpendicualar to the plane of incidence (*Hz*, *Ex*, and *Ey* are nonzero).

#### **3.2 Application of Fourier factorization**

Propagation in a 1D periodic medium, whose inhomogeneity along the *y*-axis is described by the relative permittivity function *ε*(*y*), is governed by Maxwell's equations (choosing the coordinate system uniform along the *z*-axis, i.e., *∂<sup>z</sup>* = 0)

$$
\partial\_{\mathbf{y}} E\_{\mathbf{z}} = -i\omega\mu\_{0} H\_{\mathbf{x}}, \quad \partial\_{\mathbf{x}} E\_{\mathbf{z}} = i\omega\mu\_{0} H\_{\mathbf{y}}, \quad \partial\_{\mathbf{x}} H\_{\mathbf{y}} - \partial\_{\mathbf{y}} H\_{\mathbf{x}} = i\omega\varepsilon\_{0} \varepsilon\_{(\mathcal{L})} E\_{\mathbf{z}}.\tag{31}
$$

$$
\partial\_{\mathcal{Y}} H\_{\mathcal{Z}} = i\omega \varepsilon\_0 \varepsilon\_{(\mathcal{L})} E\_{\mathcal{X}}, \quad \partial\_{\mathcal{X}} H\_{\mathcal{Z}} = -i\omega \varepsilon\_0 \varepsilon\_{(\mathcal{I})} E\_{\mathcal{Y}}, \quad \partial\_{\mathcal{X}} E\_{\mathcal{Y}} - \partial\_{\mathcal{Y}} E\_{\mathcal{X}} = -i\omega \mu\_0 H\_{\mathcal{Z}} \tag{32}
$$

where the first and the second row corresponds to the TE and the TM polarization, respectively. The formal labels of the relative permittivity, (L) and (I), means that the multiplication with the corresponding component of the electric field will be (according to the factorization rules) treated by either the Laurent rule (L) or the inverse rule (I). This is due to the fact that *Ez* and *Ex* are continuous functions (because they are tangential to the discontinuities of *ε*), whereas *Ey* is the component normal to the discontinuities and hence discontinuous (while the product *εEz* is continuous). By separating one field in each case, we obtain the wave equations

$$\begin{aligned} \begin{array}{cccc} \begin{array}{cccc} \begin{array}{c} \begin{array}{c} \begin{array}{c} \begin{array}{c} \end{array} \end{array} \end{array} \end{array} \begin{array}{c} \begin{array}{c} \begin{array}{c} \begin{array}{c} \begin{array}{c} \end{array} \end{array} \end{array} \end{array} \begin{array}{c} \begin{array}{c} \begin{array}{c} \begin{array}{c} \text{\$\boldsymbol{\sigma}\_{\mathsf{x}}^{2}}\end{array} \end{array} \begin{array}{c} \begin{array}{c} \text{\$\boldsymbol{\sigma}\_{\mathsf{x}}^{2}\$}\end{array} \end{array} \begin{array}{c} \begin{array}{c} \text{\$\boldsymbol{\tau}\_{\mathsf{x}}^{2}}\end{array} \end{array} \end{array} \begin{array}{c} \begin{array}{c} \begin{array}{c} \text{\$\boldsymbol{\tau}\_{\mathsf{x}}^{2}}\end{array} \end{array} \end{,} \begin{array}{c} \begin{array}{c} \begin{array}{c} \text{\$\boldsymbol{\tau}\_{\mathsf{x}}^{2}}\end{array} \end{array} \end{,} \begin{array}{c} \begin{array}{c} \begin{array}{c} \text{\$\boldsymbol{\tau}\_{\mathsf{x}}^{2}}\end{array} \end{array} \end{,} \end{,} \end{aligned} \end{aligned}$$

or, after expanding the permittivity or impermittivity and the fields into the Fourier series, the matrix formulae

$$\left[\mathbb{E}\right]^{-1}\left(p\_0^2 + \mathbf{q}^2\right)\left[E\_z\right] = \frac{\omega^2}{c^2}\left[E\_z\right], \quad \text{(TE)}\tag{35}$$

$$\mathbf{q}\left(\left[\frac{1}{\varepsilon}\right]p\_0^2 + \mathbf{q}\left[\varepsilon\right]^{-1}\mathbf{q}\right)[H\_{\overline{z}}] = \frac{\omega^2}{c^2}[H\_{\overline{z}}].\tag{36}$$

#### **3.3 Consistency with Yeh's theory for the small-period limit**

For the small-period limit the validity of these rules can be analytically verified by treating the periodic structure as alternation of two homogeneous layers, where we use the boundary conditions for the continuity of the tangential electric and magnetic fields on all interfaces. Now the function *ε* is assumed constant (*ε<sup>I</sup>* or *εII* ) within each layer of the thickness *d*<sup>1</sup> or *d*2. According to the geometry in Fig. 1, the field in the *j*th layer has the dependence

$$E\_z^{(j)}(x,y) = e^{-ik\_0 p\_0 x} [a\_j e^{-ik\_0 q\_j (y-y\_j)} + b\_j e^{ik\_0 q\_j (y-y\_j)}],\tag{37}$$

where *q<sup>j</sup>* =(*ε<sup>j</sup>* − *p* 2 0 ) 1/2 , with *ε<sup>j</sup>* = *ε<sup>I</sup>* for odd *j* and *ε<sup>j</sup>* = *εII* for even *j*. It is coupled with the field in the next layer by the matrix equation

$$
\begin{bmatrix} a\_{j+1} \\ b\_{j+1} \end{bmatrix} = \mathbf{C}\_{j} \begin{bmatrix} P\_{j} & 0 \\ 0 & \frac{1}{P\_{j}} \end{bmatrix} \begin{bmatrix} a\_{j} \\ b\_{j} \end{bmatrix}, \quad \mathbf{C}\_{j} = \begin{bmatrix} a\_{j} \ \beta\_{j} \\ \beta\_{j} \ a\_{j} \end{bmatrix}, \tag{38}
$$

where *α<sup>j</sup>* = 1 2 (1 + *qj*/*qj*+<sup>1</sup> ), *β<sup>j</sup>* = 1 2 (1 − *qj*/*qj*+<sup>1</sup> ) for the TE polarization, or *α<sup>j</sup>* = 1 2 (1 + *εjqj*+1/*εj*+1*q<sup>j</sup>* ), *β<sup>j</sup>* = 1 2 (1 − *εjqj*+1/*εj*+1*q<sup>j</sup>* ) for the TM polarization, and *P<sup>j</sup>* = *e* −*ik*0*qjd<sup>j</sup>* for both polarizations. Obviously, *q<sup>j</sup>* =(*ε<sup>j</sup>* − *p* 2 0 ) 1/2 , assuming the *e* −*ik*<sup>0</sup> *p*<sup>0</sup> *x* factor of all fields.

Applying the small-period approximation (*P<sup>j</sup>* ) <sup>±</sup><sup>1</sup> <sup>≈</sup> <sup>1</sup> <sup>∓</sup> *ik*0*qjd<sup>j</sup>* to the problem of propagation through the whole period,

$$
\begin{bmatrix} a\_3 \\ b\_3 \end{bmatrix} = \boldsymbol{\Omega} \begin{bmatrix} a\_1 \\ b\_1 \end{bmatrix}, \quad \boldsymbol{\Omega} = \mathbf{C}\_2 \begin{bmatrix} P\_2 & 0 \\ 0 & \frac{1}{\overline{P\_2}} \end{bmatrix} \mathbf{C}\_1 \begin{bmatrix} P\_1 & 0 \\ 0 & \frac{1}{\overline{P\_1}} \end{bmatrix}, \tag{39}
$$

yields the eigenvalues of the **Ω** operator

$$
\Omega\_{\pm} = 1 \pm ik\_0 \Lambda \sqrt{\varepsilon\_0 - p\_{0'}^2} \quad \text{(TE)} \tag{40}
$$

$$
\Omega\_{\pm} = 1 \pm ik\_0 \Lambda \sqrt{(1 - p\_0^2 \eta\_0) \varepsilon\_0}. \quad \text{(TM)} \tag{41}
$$

*ε*<sup>0</sup> = *wε<sup>I</sup>* +(1 − *w*)*εII* , (42)

$$
\eta\_0 = w \eta\_I + (1 - w) \eta\_{II} \tag{43}
$$

Assuming a Floquet mode with the eigenvalues of **<sup>Ω</sup>** being <sup>Ω</sup><sup>±</sup> <sup>=</sup> *<sup>e</sup>* <sup>±</sup>*ik*0*q*0<sup>Λ</sup> <sup>≈</sup> <sup>1</sup> <sup>±</sup> *ik*0*q*0Λ,we see that in the small-period limit the periodic structure behaves as a homogeneous anisotropic medium with the *y*-component of the normalized wave vector *q*<sup>0</sup> =(*ε*<sup>0</sup> − *p* 2 0 ) 1/2 for the TE polarization and *q*<sup>0</sup> =[(1 − *p* 2 0 *η*0)*ε*0] 1/2 for the TM polarization. These formulae are identical with Equations 35 and 36 if we retain only the 0th element of all matrices. This is a very interesting disclosure that the results obtained by Yeh and coauthors already in 1970s are consistent with the extensive, more general research carried out in 1990.

## **3.4 Numerical example**

Example of the comparison of applying the correct Fourier factorization rules with applying the opposite ones is shown in Fig. 2 for both polarizations. The normalized eigenfrequency *ω*Λ/2*πc* of the first band is displayed according to *N*; the structure is made as two alternating layers of the equal thicknesses 500 nm (Λ = 1000 nm) and permittivities *ε<sup>I</sup>* = 3 and *εII* = 1. The wave vector is chosen k =(0.5*π*/Λ)[1, 1, 0].

Fig. 2. Convergence properties of the correct factorization compared with the opposite one

# **4. Two-dimensional photonic crystals with circular elements**

#### **4.1 Geometry of the problem**

Fig. 3 displays the geometrical arrangement of a 2D photonic crystal made as bi-periodic alternation of rods or holes with a cylindrical cross-section. Instead of a single vector of periodicity we now have the plane of periodicity (determined by two vectors of periodicities defining a unit cell), which here coincides with the *xy* plane. For simplicity we choose the incidence direction in the plane of periodicity, so that we can again distinguish between two independent polarizations. The plane of incidence is now determined along the propagation wave vector k = *k*0[*p*0, *q*0,0] and perpendicular to the plane of periodicity. The TE polarization has now H along the *z* axis, which is now the more difficult case (with nonzero

with

Fig. 3. Geometry of a 2D photonic crystal made as bi-periodic alternation of rods or holes

*Hz*, *Ex*, *Ey*), while the TM polarization has E along the *z* axis (nonzero *Ez*, *Hx*, *Hy*), which is now the simple case.

Fig. 4. Two examples of 2D periodic arrangements; below are the first Brillouin zones in the reciprocal space

As an example, in Fig. 3 we have chosen the *xz* plane as the plane of incidence, but we could choose any plane parallel with the *z* axis provided that we want to treat the TE and TM polarizations independently. In general, the plane of incidence is determined by the k vector in the reciprocal space as displayed in Fig. 4, whose particular symmetry points are denoted Γ, X, M, or K according to the corresponding periodicity.

In this section we study a 2D photonic crystal composed of infinite cylinders with a circular cross-section with either square [Fig. 4(a)] or hexagonal [Fig. 4(b)] periodicity. For the square symmetry the unit cell has the dimensions Λ*<sup>x</sup>* = Λ*<sup>y</sup>* = Λ, and for the hexagonal symmetry

it can be chosen as an Λ-by-Λ √ 3 rectangle. The corresponding first Brillouin zones of the reciprocal space are depicted in the bottom part of Fig. 4.

Maxwel's equation are again

$$
\partial\_y E\_z = -i\omega\mu\_0 H\_\Sigma, \quad \partial\_\mathbf{x} E\_\mathbf{z} = i\omega\mu\_0 H\_\mathbf{y}, \quad \partial\_\mathbf{x} H\_\mathbf{y} - \partial\_\mathbf{y} H\_\mathbf{x} = i\omega\varepsilon\_0 \varepsilon E\_\mathbf{z} \tag{44}
$$

$$
\partial\_y H\_2 = i\omega \varepsilon\_0 \varepsilon E\_\chi \quad \partial\_x H\_2 = -i\omega \varepsilon\_0 \varepsilon E\_y \quad \partial\_x E\_y - \partial\_y E\_x = -i\omega \mu\_0 H\_2. \tag{45}
$$

However, now we cannot put there the labels (L) or (I) for the factorization type as easily as above, because the discontinuities of the permittivity function are now mixed among the components of E.

Assuming a hypothetical anisotropy of the relative permittivity function, we define a scaled electrical displacement D ˜ ,

$$
\begin{bmatrix} \tilde{D}\_{\mathcal{X}} \\ \tilde{D}\_{\mathcal{Y}} \end{bmatrix} = \mathfrak{e} \begin{bmatrix} E\_{\mathcal{X}} \\ E\_{\mathcal{Y}} \end{bmatrix} = \begin{bmatrix} \varepsilon\_{\mathcal{X}\mathcal{X}} \ \varepsilon\_{\mathcal{Y}\mathcal{Y}} \\ \varepsilon\_{\mathcal{Y}\mathcal{X}} \ \varepsilon\_{\mathcal{Y}\mathcal{Y}} \end{bmatrix} \begin{bmatrix} E\_{\mathcal{X}} \\ E\_{\mathcal{Y}} \end{bmatrix}, \tag{16}
$$

$$
\mathbf{D}\_2 = \varepsilon\_{22} E\_2 = \varepsilon\_{(\mathcal{L})} E\_2 \quad \text{ (TM)}\tag{47}
$$

where *εzz* is obviously the only component of the permittivity tensor for which we can use the Laurent rule.

Defining also a 2-by-2 matrix of electrical impermittivity η = ε <sup>−</sup><sup>1</sup> helps in the formulation of the wave equations

$$(-\partial\_{\mathbf{y}}\eta\_{\mathbf{xx}}\partial\_{\mathbf{y}} + \partial\_{\mathbf{x}}\eta\_{\mathbf{yx}}\partial\_{\mathbf{y}} + \partial\_{\mathbf{y}}\eta\_{\mathbf{xy}}\partial\_{\mathbf{x}} - \partial\_{\mathbf{x}}\eta\_{\mathbf{yy}}\partial\_{\mathbf{x}})H\_{\mathbf{z}} = \frac{\omega^2}{c^2}H\_{\mathbf{z}}, \quad \text{(TE)}\tag{48}$$

$$-\left(\partial\_x^2 + \partial\_y^2\right) E\_2 = \frac{\omega^2}{c^2} \varepsilon\_{(\mathcal{L})} E\_{2,\prime} \quad \text{(TM)}\tag{49}$$

where *ηjk* are the components of the electrical impermittivity. For the simplicity of the TM polarization case we below focus our attention only to the TE polarization.

#### **4.2 Methods of Fourier factorization**

In this section we compare several models corresponding to different factorization approaches.

#### **4.2.1 Elementary (Cartesian) method (Model A)**

First, Model A assumes the solution in the basis of the ˆx and ˆy polarizations uniform within the periodic cell, where in accordance with Ho et al. (1990) we choose the Laurent rules

$$\begin{bmatrix} \mathbf{\tilde{D}\_{\mathbf{x}}} \end{bmatrix} = \begin{bmatrix} \varepsilon \mathbf{E}\_{\mathbf{x}} \end{bmatrix} = \begin{bmatrix} \mathbf{\tilde{c}} \end{bmatrix} \begin{bmatrix} \mathbf{E}\_{\mathbf{x}} \end{bmatrix} \tag{50}$$

$$[\tilde{D}\_y] = [\varepsilon E\_y] = [\mathbb{[\varepsilon]}][E\_y]. \tag{51}$$

The components of the electric impermittivity in Equation 48 then becomes [[*ε*]] −1 for the cases of *ηxx*, *ηyy*, and zero for the cases of *ηxy*, *ηyx*,or

$$\begin{array}{c} \llbracket \boldsymbol{\eta} \rrbracket \rrbracket\_{\mathcal{A}} = \begin{bmatrix} \llbracket \boldsymbol{\varepsilon} \rrbracket \rrbracket^{-1} \begin{bmatrix} \llbracket \boldsymbol{0} \rrbracket \\\llbracket \boldsymbol{\varepsilon} \rrbracket \rrbracket^{-1} \end{bmatrix}. \end{array} . \tag{52}$$

For illustration we show the distribution of the first basis polarization vector (identical with the constant vector ˆx) in Fig. 5(a), where the black circle denotes the element boundary (the permittivity discontinuity).

**4.2.2 Normal vector method (Model B)**

Fig. 5. Distribution of the basis polarization vector u for the factorization models

According to the factorization theorems, neither the Laurent rule nor the inverse rule is correct for both products in Equations 50 and 51, because both pairs of functions have concurrent discontinuities and both products *D* ˜ *x* and *D* ˜ *y* are discontinuous as well. On the other hand, by an appropriate change of the polarization bases at all points (using a space-dependent Jones matrix transform **F**),

$$
\begin{bmatrix} E\_x \\ E\_y \end{bmatrix} = \mathbf{F} \begin{bmatrix} E\_{\mathcal{U}} \\ E\_{\mathcal{U}} \end{bmatrix} \tag{53}
$$

we can treat independently the normal (*u*) and tangential (*v*) components of the fields by the correct rules,

$$[\tilde{D}\_{\mu}] = [\![1/\varepsilon]\!]^{-1}[E\_{\mu}]\_{\prime} \tag{54}$$

$$[\![\![\!D\_{\upsilon}]\!] = [\![\varepsilon]\!][E\_{\upsilon}].\tag{55}$$

The field components *Eu*, *D* ˜ *u* are normal to the discontinuities of the relative permittivity function, while *Ev*, *D* ˜ *v* are tangential. The factorization rules used in Equations 54 and 55 are justified simply because *Ev* and *D* ˜ *u* are continuous.

A suitable distribution of the matrix **F** within the periodic cell can obviously be the rotation

$$\mathbf{F} = \begin{bmatrix} \cos \phi - \sin \phi \\ \sin \phi & \cos \phi \end{bmatrix}'\tag{56}$$

where the polar angle *φ*(*x*, *y*) is in the first cell distributed according to the polar coordinates *rei<sup>φ</sup>* = *x* + *iy*, and then periodically repeated over the entire 2D space. This enables defining the matrices [[*c*]], [[*s*]] from the corresponding 2D-periodic functions *c* = cos *φ*, *s* = sin *φ*.

Let u and v be the two columns of the matrix **F**, both being mutually orthogonal basis vectors of linear polarization. From the above definitions we see that u is a polarization vector normal to the structure discontinuities, whereas v is tangential. In Fig. 5(b) we show the distribution of u within the periodic cell. The basis polarization vectors are constant along the lines of the constant azimuth (*φ* = const) and rotate as *φ* increases. It is obvious that the matrix function **F**(*x*, *y*) has no discontinuities concurrent with the electric field, so that we can use both Laurent and inverse rules for the transformation of polarization, e.g.,

$$
\begin{bmatrix} \begin{bmatrix} E\_{\mathcal{X}} \end{bmatrix} \\ \begin{bmatrix} E\_{\mathcal{Y}} \end{bmatrix} \end{bmatrix} = \begin{bmatrix} \mathbf{F} \end{bmatrix} \begin{bmatrix} \begin{bmatrix} E\_{\mathcal{U}} \end{bmatrix} \\ \begin{bmatrix} E\_{\mathcal{U}} \end{bmatrix} \end{bmatrix} \tag{57}
$$

$$\begin{array}{c} \begin{bmatrix} \mathbf{F} \end{bmatrix} = \begin{bmatrix} \begin{bmatrix} \mathbf{c} \end{bmatrix} \begin{bmatrix} -s \end{bmatrix} \\\ \begin{bmatrix} \mathbf{s} \end{bmatrix} \begin{bmatrix} \mathbf{c} \end{bmatrix} \end{array} \end{array} \tag{58}$$

Combining Equations 54, 55, and 57 yields

$$
\begin{bmatrix} \begin{bmatrix} E\_X \\ E\_y \end{bmatrix} \end{bmatrix} = \begin{bmatrix} \mathbf{F} \end{bmatrix} \begin{bmatrix} \begin{bmatrix} \frac{1}{\varepsilon} \end{bmatrix} \begin{bmatrix} \mathbf{0} \end{bmatrix} \\ \begin{bmatrix} \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{0} \end{bmatrix}^{-1} \begin{bmatrix} \mathbf{F}^{-1} \end{bmatrix} \begin{bmatrix} \begin{bmatrix} \tilde{D}\_x \\ D\_y \end{bmatrix} \end{bmatrix} \end{bmatrix} \tag{59}
$$

from where we derive the electric impermittivity in the reciprocal space (corresponding to Model B)

$$\begin{aligned} \begin{bmatrix} \mathfrak{h} \end{bmatrix}\_{\mathbb{B}} &= \begin{bmatrix} \mathbb{F} \end{bmatrix} \begin{bmatrix} \begin{bmatrix} \frac{1}{\varepsilon} \\ \mathbb{0} \end{bmatrix} \begin{bmatrix} \mathbb{O} \\ \mathbb{e} \end{bmatrix}^{-1} \end{bmatrix} \begin{bmatrix} \mathbb{F}^{-1} \\ \end{bmatrix} \\ &= \begin{bmatrix} \mathbb{I} \begin{bmatrix} \frac{1}{\varepsilon} \\ \mathbb{I} \end{bmatrix} \begin{bmatrix} \varepsilon^{2} \\ \end{bmatrix} + \begin{bmatrix} \mathbb{I} \end{bmatrix} \begin{bmatrix} \mathbb{I}^{2} \end{bmatrix} \begin{bmatrix} \mathbb{I} \end{bmatrix} \begin{bmatrix} \mathbb{I} \end{bmatrix} \begin{bmatrix} \mathbb{C} \\ \mathbb{S} \end{bmatrix} \begin{bmatrix} \mathbb{C} \\ \mathbb{S} \end{bmatrix} - \begin{bmatrix} \mathbb{I} \end{bmatrix}^{-1} \begin{bmatrix} \mathbb{C} \\ \mathbb{S} \end{bmatrix} \\ &= \begin{bmatrix} \mathbb{I} \end{bmatrix} \begin{bmatrix} \mathbb{I} \mathbb{S} \\ \mathbb{I} \end{bmatrix} \begin{bmatrix} \mathbb{C} \\ \mathbb{S} \end{bmatrix} - \begin{bmatrix} \mathbb{I} \end{bmatrix} \begin{bmatrix} \mathbb{I} \end{bmatrix} \begin{bmatrix} \mathbb{S} \\ \mathbb{S} \end{bmatrix} \begin{bmatrix} \mathbb{S}^{2} \\ \mathbb{S} \end{bmatrix} + \begin{bmatrix} \mathbb{I} \end{bmatrix}^{-1} \begin{bmatrix} \mathbb{C} \\ \mathbb{C} \end{bmatrix} \end{bmatrix} \end{aligned} \tag{60}$$

whose components are immediately applicable to Equation 48.

### **4.2.3 Method with elliptical polarization bases (Model C)**

The above approach (Model B) only deals with linear polarizations and thus suffers from the fact that the matrix function **F**(*x*, *y*) has a singularity at the central point of the periodic cell and other discontinuities along the cell boundaries. This slows down the convergence of the numerical implementation, as will be evidenced below. On the other hand, we can make **F** continuous by using complex functions *ξ* and *ζ* or, in other words, by defining u, v as complex vectors corresponding to generally elliptic polarizations,

$$u = \begin{bmatrix} \mathfrak{J} \\ \mathfrak{J} \end{bmatrix}, \quad v = \begin{bmatrix} -\mathfrak{J}^\* \\ \mathfrak{J}^\* \end{bmatrix} \tag{61}$$

(still orthogonal). By means of rotation *θ* and ellipticity *ǫ* we define the first basis vector

$$u = e^{i\theta} \begin{bmatrix} \cos\theta - \sin\theta \\ \sin\theta \end{bmatrix} \begin{bmatrix} \cos\epsilon \\ i\sin\epsilon \end{bmatrix} \tag{62}$$

where

$$\begin{array}{c} \square \quad \theta(r,\phi) = \phi, \\\square \quad \square \quad \begin{array}{c} \theta(r,\phi) = \phi, \\\varepsilon(r,\phi) = \begin{cases} \frac{\pi}{8} \left(1 + \cos\frac{\pi r}{R}\right) & (r \le R) \\\ \frac{\pi}{8} \left(1 + \cos\frac{\pi \left[r + D(\phi) - 2R\right]}{D(\phi) - R}\right) & (r > R) \end{cases} & \xrightarrow{\square} \end{array} \end{array} \tag{63}$$

8 1 + cos *D*(*φ*)−*R* Here *R* denotes the radius of the circular element and

$$D(\phi) = \frac{\Lambda/2}{\max(|\cos \phi|\_{\prime} | \sin \phi|)}\tag{65}$$

is the distance from the cell's center to its edge. In Equation 62 the Jones vector on the right represents a polarization ellipse (with ellipticity *ǫ*) oriented along the *x* coordinate, the matrix in the middle rotates this polarization by the azimuth *θ*, and the factor *e <sup>i</sup><sup>θ</sup>* preserves the continuity of the phase at the center and along the boundaries of the cell. This continuity can be easily checked by evaluating the limits

$$\lim\_{r \to 0} u = \lim\_{r \to D(\phi)} u = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ i \end{bmatrix} \, , \tag{66}$$

which is the vector of left circular polarization (independent of *φ*).

The distribution of the basis polarization vector u within the periodic cell is shown in Fig. 5(c). Here the azimuth of the polarization ellipse is constant along the lines coming from the cell's center, which is similar to Model B. However, the ellipticity is now zero (corresponding to linear polarization) only on the boundaries of the circular element, has the maximum value (*π*/4 for circular polarization) at the cell's center and along its boundaries, and continuously varies (with a smooth sine dependence) in the intermediate ranges. Thus we obtain a smooth and completely continuous matrix function **F**(*x*, *y*), which is analogously used to calculate the impermittivity in the reciprocal space

[[η]]<sup>C</sup> = [[ 1 *ε* ]] [[*ξξ* ∗ ]]+[[*ε*]] −1 [[*ζζ* ∗ ]], [[ 1 *ε* ]] [[*ξζ* ∗ ]] − [[*ε*]] −1 [[*ξζ* ∗ ]] [[ 1 *ε* ]] [[*ξ* ∗ *ζ*]] − [[*ε*]] −1 [[*ξ* ∗ *ζ*]], [[ 1 *ε* ]] [[*ζζ* ∗ ]]+[[*ε*]] −1 [[*ξξ* ∗ ]] . (67)

In the case of the hexagonal periodicity we define u and the other periodic quantities inside one hexagon (half the area of the rectangular unit cell) where we can use formally the same equations as above, except for

$$D(\phi) = \frac{\Lambda/2}{\max\_{n=0,\ldots,5} \left[ \cos \left( \phi - \frac{n\pi}{3} \right) \right]},\tag{68}$$

which is now the distance from the hexagon's center to its edge. Here Λ is the hexagon's shortest width (equal to the width Λ*<sup>x</sup>* of the rectangular cell).

#### **4.2.4 Modified method for densely arranged elements (Model C')**

To analyze a more complicated situation, we consider a photonic crystal with square periodicity where circular elements are densely arranged near each other, i.e., where the radius *R* is almost the half width Λ/2 of the periodic cell. Then the convergence properties of **F** becomes worse, which affects all the derived quantities. For this reason we again redefine the polarization distribution. For the modified Model C' we define u to be still same inside the circle (*r* < *R*), but different outside. Assuming the rotation and ellipticity along the boundary of the square cell

$$\begin{array}{ccccc}\hline\\\hline\\\hline\end{array}\quad\begin{array}{ccccc}\hline\theta\_{\mathsf{b}}(\phi)&=\theta\_{\mathsf{b}}(D(\phi),\phi)&=\frac{\pi}{2}\texttt{round}\left(\phi/\frac{\pi}{2}\right),&\boxed{\qquad\smile\,\,\,}\\\hline\\\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{array}\quad\begin{array}{ccccc}\hline\hline\end{$$

$$\epsilon\_{\rm b}(\phi) = \overline{\epsilon \left( D(\phi), \phi \right)} = \frac{\pi}{8} (1 - \cos 4\phi) \tag{70}$$

(where "**round**" denotes rounding towards the nearest integer), we define the rotation and ellipticity outside the circle (*r* > *R*)as

$$\theta(r,\phi) = \frac{1}{2}\left\{\theta\_{\mathsf{b}}(\phi) + \phi + [\theta\_{\mathsf{b}}(\phi) - \phi]\cos\frac{\pi[r + D(\phi) - 2\mathsf{R}]}{D(\phi) - \mathsf{R}}\right\},\tag{71}$$

$$\varepsilon(r,\phi) = \frac{\varepsilon\_{\rm b}(\phi)}{2} \left\{ 1 + \cos \frac{\pi [r + D(\phi) - 2R]}{D(\phi) - R} \right\}. \tag{72}$$

Assuming otherwise the same Equations 59, 61, 62, and 65, we obtain for [[η]]C′ formally the same matrix as in Equation 67, except that the functions *ξ* and *ζ* are now derived from different azimuth and ellipticity distributions of u. Note that u is again continuous along the cell's boundaries; to evaluate its precise limits [when *x* →±*D*(0), *y* = const or *y* →±*D*(*φ*/2), *x* = const] would now be more complicated. The distribution of the basis polarization vector u within the periodic cell is depicted in Fig. 5(d) together with dimensions.

#### **4.3 Numerical examples**

We examine the numerical performances of all factorization models presented in Section 4.2 on three samples of 2D photonic crystals, for which we calculate the eigenfrequencies *ωκ* (where the band number *κ* = 1 stands for the lowest eigenfrequency, *κ* = 2 for the second lowest, etc.) and the corresponding eigenvectors [*Hz*]*<sup>κ</sup>* of Equation 48. All convergences will be presented according to the maximum Fourier harmonics retained inside the periodic medium, which will be kept same for the *x* and *y* directions (*M* = *N*).

First, Sample S1 is a square array of cylindrical rods of the circular cross-section with the diameter 2*R* = 500 nm, square period Λ = 1000 nm, relative permittivity of the rods *ε*<sup>1</sup> = 9, and relative permittivity of the surrounding medium corresponding to vacuum (*ε*<sup>2</sup> = 1). Its dispersion relation is displayed in Fig. 6(a). Similarly, for Sample S2 we assume exactly the same parameters except the diameter of the rods, now being 2*R* = 900 nm. This corresponds to densely arranged elements (the distance between two adjacent rods is only 100 nm). Finally, for Sample H we consider a hexagonal array of cylindrical holes of the circular cross-section with the diameter 2*R* = 600 nm, hexagonal periodicity Λ = 1000 nm (corresponding to the rectangular cell of the dimensions Λ*<sup>x</sup>* = 1 *μ*m, Λ*<sup>y</sup>* = √ 3 *μ*m), relative permittivity of the holes corresponding to vacuum (*ε*<sup>1</sup> = 1), and relative permittivity of the substrate medium (surrounding holes) *ε*<sup>2</sup> = 12. The dispersion relation of Sample H is displayed in Fig. 6(b).

For our analysis we choose the eigenmodes Γ2 and Γ3 of Sample S1, the eigenmode X3 of Sample S2, and the eigenmode M1 of Sample H, where the letter denotes a point of symmetry

Fig. 6. Photonic band structures of two samples

Fig. 7. Amplitude distributions of the scaled magnetic field |*Hz*| within the unit cell of four chosen eigenmodes

and the number corresponds to the band. The amplitude distributions of the modes are displayed in Fig. 7, and the corresponding convergence properties of the eigenfrequencies calculated by the above described models are shown in Fig. 8 (Model C' is only compared for Sample S2).

The result of a careful comparison of the numerical efficiencies of all the factorization models can be summarized as follows. Models B, C, and C' always converge considerably faster than Model A. Model C converges faster (having usually one order higher precision) than Model B with two exceptions. The first exception is the case such as in Fig. 7 (S1-Γ2) and Fig. 8(a), where the discontinuities of the polarization transformation matrix **F** coincide with a nearly zero amplitude of the field, so that the discontinuities do not manifest themselves. The second exception is the case such as Fig. 7 (S2-X3) and Fig. 8(d), where the elements are densely arranged which causes rapid variations of the ellipticity between two adjacent elements (which are very close to each other); this requires more Fourier components than the weak discontinuity of the linear polarization u in Model B. The problem is solved in Model C', which obviously converges fastest among all the four models applied to Sample S2.

Fig. 8. Convergences of the Fourier factorization methods for three samples

# **5. Two-dimensional photonic crystals with non-circular elements**

In this section we briefly show numerical efficiencies of three factorization models derived above (Models A, B, and C) applied to 2D photonic crystals made as long elements of other shapes, namely rods with the square cross-section and tubes with the ring and split-ring cross-sections, arranged with the square periodicity. For all the three samples we choose the permittivity of the elements *ε* = 9 and permittivity of vacuum *ε* = 1 for the surrounding medium.

# **5.1 Periodic rods with the square cross-section**

For the photonic crystals made of square-sectioned rods, the period is chosen Λ = 1000 nm, and the width of the square *d* = 600 nm. The distribution of the basis polarization vector u, analogously to Section 4, is displayed in Fig. 9 for all the three compared models.

For Model A, of course, u = x ˆ as visible in Fig. 9(a). For Model B, we divide the unit cell into four parts by its diagonals (the lines *x* = *y* and *x* = −*y*). The vector u(*φ*), depending only

Fig. 9. Polarization distribution of the factorization methods for periodic squares

on the azimuthal angle, is set ˆ<sup>x</sup> for *<sup>φ</sup>* ∈−*<sup>π</sup>* 4 , *π* 4 ) ∪ 3*π* 4 , 5*π* 4 ), and ˆy for the remaining angles. This means that u is always linear, perpendicular to the permittivity discontinuities, and its discontinuities (along the lines *x* = *y* and *x* = −*y*) have no concurrent discontinuities with the permittivity function except those at the four corners of the square. Hence, Model B here fulfills nearly the same conditions for the application of the factorization rules as demanded in Section 4 for circular elements.

For Model C we divide the unit cell into four areas in the same manner. This time, however, the basis vector u(*r*, *φ*) depends on both polar coordinates and the corresponding polarization is in general elliptic. In analogy with Section 4 we want u to be perpendicular to the permittivity discontinuities and to remove its discontinuities as much as possible. The most simple way how to do this, although the discontinuities will not be completely removed, is to set the azimuth *<sup>θ</sup>*(*φ*) of the polarization to zero for *<sup>φ</sup>* ∈−*<sup>π</sup>* 4 , *π* 4 ) ∪ 3*π* 4 , 5*π* 4 ), and *π* 2 for the remaining angles, and to use Equation 64 for the ellipticity distribution, where *D*(*φ*) now corresponds to the square element.

Fig. 10. Convergence of the factorization methods for rods with the square cross-section

The numerical efficiencies of all the three models are displayed in Fig. 10 for the 3rd band of the symmetry point Γ. As clearly visible, Models B and C converge considerably faster than Model A, with Model C being a little better. Although the improvement from Model B towards Model C is not so distinct, it could be improved by a more careful choice of a completely continuous distribution of the polarization basis.

## **5.2 Periodic hollow cylinders**

For the photonic crystals made of periodic hollow cylinders (tubes with the symmetric ring cross-section) we choose the period Λ = 1000 nm, the inner diameter (the diameter of the inner circular hole) *R*<sup>1</sup> = 400 nm, and the outer diameter *R*<sup>2</sup> = 680 nm. The distribution of the basis polarization vector u for all the three compared models is displayed in Fig. 11.

Fig. 11. Polarization distribution of the factorization methods for hollow cylinders

For Model A again u = x ˆ. For both Models B and C we use the same polarization distributions as in Section 4 for the inner area (*r* < *R*<sup>1</sup> ) and for the outer area (*r* > *R*2). For the annulus area (*R*<sup>1</sup> < *r* < *R*2) we simply choose the polarization distribution same as that in Model B. These distributions are obviously the most straightforward analogies of the distributions used in Section 4 for circular elements.

The numerical efficiencies of all the three models are displayed in Fig. 12 for the 3rd band of the symmetry point Γ. As clearly visible, Models B and C converge considerably faster than Model A, but now Model C exhibits no improvement against Model B. This is because the discontinuities of u at the center and along the boundaries of the periodic cell quite well coincide with the zero amplitude of the mode, as visible in the inset of Fig. 12, so that the discontinuities do not manifest themselves in the calculations.

# **5.3 Periodic split hollow cylinders**

For the photonic crystals made of split hollow cylinders (tubes with an asymmetric, split ring cross-section) we choose the period Λ = 1000 nm, the inner diameter (the diameter of the inner semi-circular hole) *R*<sup>1</sup> = 600 nm, the outer diameter *R*<sup>2</sup> = 720 nm, and the relative azimuthal length of the ring *w<sup>φ</sup>* = 0.9 (where *w<sup>φ</sup>* = 1 means the complete, symmetric ring). The distributions of the basis polarization vector u for all the three compared models are displayed in Fig. 13.

Fig. 12. Convergence of the factorization methods for hollow cylinders

Fig. 13. Polarization distribution of the factorization methods for split hollow cylinders

Since the area of splitting is quite small compared to the full area of the unit cell, we have chosen the polarization distributions of all the three models exactly same as for the symmetric rings in Section 5.2. The condition for normal u and tangential v is not satisfied on the surface proportional to the length 2 × 120 = 240 nm, but is satisfied on the surface proportional to 0.9 × 2*π*(*R*<sup>1</sup> + *R*2) ≈ 7 464 nm, which is 30 times higher area, justifying this negligence. Of course, for modes with most of electromagnetic field resonating in the critical area this approximation would be insufficient.

The numerical efficiencies of all the three models are displayed in Fig. 14, again for the 3rd band of the symmetry point Γ. We can describe these performances by exactly the same conclusion as for the symmetric rings in Section 5.2. Models B and C converge similarly and both considerably faster than Model A. The discontinuities of u at the center and along the boundaries of the periodic cell again coincide with the zero amplitude of the mode, though with some deviations visible in the inset of Fig. 14.

Fig. 14. Convergence of the factorization methods for split hollow cylinders

# **6. Anisotropic photonic crystals**

In this section we will briefly demonstrate the application of the factorization rules to 2D photonic crystals made of anisotropic materials, again with the plane of incidence parallel to the *z* axis (with geometry of Fig. 3). Unlike the isotropic crystals, now the TE and TM polarizations are not separable. Instead of the scalar permittivity we define and expand the components of the relative permittivity tensor function

$$\varepsilon\_{jk}(x,y) = \sum\_{m,n=-\infty}^{+\infty} \varepsilon\_{jk,mn} e^{-i(mpx+nqy)}\,,\tag{73}$$

where *εjk*,*mn* are the Fourier coefficients. The wave equation for a generally anisotropic medium, now described by the permittivity tensor

$$\mathbf{c} = \begin{bmatrix} \varepsilon\_{xx} \ \varepsilon\_{xy} \ \varepsilon\_{xz} \\ \varepsilon\_{yx} \ \varepsilon\_{yy} \ \varepsilon\_{yz} \end{bmatrix}. \tag{74}$$
  $\mathbf{c} \text{ is the operator equation}$  
$$\mathbf{c} \text{ is the operator equation}$$
 
$$\mathbf{\hat{C}E} = \frac{\omega^2}{c^2} \mathbf{E}, \quad \mathbf{\hat{C}} = \varepsilon^{-1} \begin{bmatrix} -\mathbf{\hat{o}}\_y^2 & \mathbf{\hat{o}}\_x \mathbf{\hat{o}}\_y & \mathbf{0} \\ \mathbf{\hat{o}}\_x \mathbf{\hat{o}}\_y & -\mathbf{\hat{o}}\_x^2 & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & -\mathbf{\hat{o}}\_x^2 - \mathbf{\hat{o}}\_y^2 \end{bmatrix}. \tag{75}$$

Similarly as above, Model A assumes all components of the permittivity tensor treated by the Laurent rule, i.e.,

$$\mathbb{E}\left[\tilde{D}\_{j}\right] = \sum\_{k} \mathbb{I}\varepsilon\_{jk} \mathbb{I}\left[E\_{k}\right]. \tag{76}$$

To apply the factorization correctly, we must again separate the normal and tangential components of all fields for which we can use the correct rules. Let us define a space-dependent matrix transform **F**(*x*, *y*) so that

$$
\begin{bmatrix} E\_{\mathcal{X}} \\ E\_y \\ E\_z \end{bmatrix} = \mathbf{F} \begin{bmatrix} E\_u \\ E\_v \\ E\_z \end{bmatrix}'\tag{77}
$$

where *E<sup>u</sup>* and *E<sup>v</sup>* are the normal and tangential components of the vector E to all discontinuities of the permittivity.

Analogously to Section 4, for the case of circular elements we can choose the polarization basis distribution corresponding to Models B and C as

$$\mathbf{F} = \begin{bmatrix} \cos \phi - \sin \phi & 0 \\ \sin \phi & \cos \phi & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \text{(Model B)}\tag{78}$$

$$\mathbf{F} = \begin{bmatrix} \mathfrak{J} & -\mathfrak{J}^\* & 0 \\ \mathfrak{J} & \mathfrak{J}^\* & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \text{(Model } \mathbf{C} \text{)}\tag{79}$$

where the matrix elements have the same meaning as in Section 4. In the new coordinates we can write

$$
\begin{bmatrix} \tilde{D}\_{\boldsymbol{u}} \\ \tilde{D}\_{\boldsymbol{v}} \\ \tilde{D}\_{\boldsymbol{z}} \end{bmatrix} = \begin{bmatrix} \varepsilon\_{\ell\mu\boldsymbol{u}} \ \varepsilon\_{\ell\nu\boldsymbol{v}} \ \varepsilon\_{\ell\boldsymbol{z}} \\ \varepsilon\_{\varepsilon\mu\boldsymbol{u}} \ \varepsilon\_{\nu\boldsymbol{v}} \ \varepsilon\_{\ell\boldsymbol{v}} \\ \varepsilon\_{\varepsilon\mu\boldsymbol{u}} \ \varepsilon\_{\varepsilon\nu} \ \varepsilon\_{\varepsilon\boldsymbol{z}} \end{bmatrix} \begin{bmatrix} E\_{\boldsymbol{u}} \\ E\_{\boldsymbol{v}} \\ E\_{\boldsymbol{z}} \end{bmatrix} . \tag{80}
$$

Now let us separate two sets of quantities, those which are continuous (*D* ˜ *u*, *Ev*, *Ez*) and those which are not continuous (*D* ˜ *v*, *D* ˜ *z*, *Eu*) to the discontinuities of the permittivity. Expressing the second set according to the first one yields

$$
\begin{bmatrix}
\mathbf{E}\_{u} \\
\mathbf{\tilde{D}}\_{v} \\
\mathbf{\tilde{D}}\_{z}
\end{bmatrix} = \mathbf{G} \begin{bmatrix}
\tilde{D}\_{u} \\
\mathbf{E}\_{v} \\
\mathbf{\tilde{D}}\_{z}
\end{bmatrix}, \quad \mathbf{G} = \begin{bmatrix}
\frac{\mathcal{I}\_{\rm luv}}{\mathcal{E}\_{\rm luv}} - \frac{\mathcal{E}\_{\rm luv}}{\mathcal{E}\_{\rm luv}} & -\frac{\mathcal{E}\_{\rm luv}}{\mathcal{E}\_{\rm luv}} \\
\frac{\mathcal{E}\_{\rm luv}}{\mathcal{E}\_{\rm luv}} \varepsilon\_{\rm vu} - \frac{\mathcal{E}\_{\rm luv} \mathcal{E}\_{\rm luv}}{\mathcal{E}\_{\rm luv}} \varepsilon\_{\rm v2} - \frac{\mathcal{E}\_{\rm luv} \mathcal{E}\_{\rm luv}}{\varepsilon\_{\rm vu}} \\
\end{bmatrix} \tag{81}
$$
  $\text{which we can simply use the Laurent rule}$ 

for which we can simply use the Laurent rule,

$$
\begin{bmatrix} \left[ E\_{\mathcal{U}} \right] \\ \left[ \left[ \mathcal{D}\_{\upsilon} \right] \right] \\ \left[ \left[ \mathcal{D}\_{\mathcal{Z}} \right] \right] \end{bmatrix} = \left[ \mathbf{G} \right] \begin{bmatrix} \left[ \tilde{D}\_{\mathcal{U}} \right] \\ \left[ E\_{\upsilon} \right] \\ \left[ E\_{\mathcal{Z}} \right] \end{bmatrix} . \tag{82}
$$

From this we express [*D* ˜ *j* ] according to [*E<sup>j</sup>* ],

$$
\begin{bmatrix} \left[\mathcal{D}\_{u}\right] \\ \left[\tilde{D}\_{v}\right] \\ \left[\tilde{D}\_{z}\right] \end{bmatrix} = \left[\mathbb{E}\_{\{uvz\}}\right] \mathbb{I}\_{\mathcal{B},\mathbb{C}} \begin{bmatrix} \left[E\_{u}\right] \\ \left[E\_{v}\right] \\ \left[E\_{\mathbb{Z}}\right] \end{bmatrix} \tag{83}
$$

with

[[ε{*uvz*} ]] B, C = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 *εuu* <sup>−</sup><sup>1</sup> , 1 *εuu* <sup>−</sup><sup>1</sup> *εuv εuu* , *εvu εuu* 1 *εuu* <sup>−</sup><sup>1</sup> , *εvv* − *εvuεuv εuu* + *εvu εuu* 1 *εuu* <sup>−</sup><sup>1</sup> *εuv εuu* , *εzu εuu* 1 *εuu* <sup>−</sup><sup>1</sup> , *εzv* − *εzuεuv εuu* + *εzu εuu* 1 *εuu* <sup>−</sup><sup>1</sup> *εuv εuu* , 1 *εuu* <sup>−</sup><sup>1</sup> *εuz εuu εvz* − *εvuεuz εuu* + *εvu εuu* 1 *εuu* <sup>−</sup><sup>1</sup> *εuz εuu εzz* − *εzuεuz εuu* + *εzu εuu* 1 *εuu* <sup>−</sup><sup>1</sup> *εuz εuu* ⎤ ⎥ ⎥ ⎥ ⎥ . (84)

Finally we obtain the matrix of permittivity in the Cartesian coordinates via the formula

$$\left\| \left[ \varepsilon\_{\{xyz\}} \right] \right\|\_{\mathrm{B},\mathrm{C}} = \left\| \mathrm{F} \right\| \left[ \varepsilon\_{\{\mu vz\}} \right] \right\|\_{\mathrm{B},\mathrm{C}} \left[ \mathrm{F}^{-1} \right] \mathrm{.}\tag{85}$$

where the index B or C corresponds to the chosen model.

# **7. Distribution of polarization basis for 3D structures**

In this section we will suggest how to create the polarization basis distribution for a 3D photonic crystal. Unlike the previous cases, we now have a fully vectorial wave Equation 1. The corresponding material equation D ˜ = *ε*(*x*, *y*, *z*)E must be changed to separate the normal and tangential components of the fields to the *ε* discontinuities, which are now surfaces in the 3D space. For simplicity we assume a 3D photonic crystal made as spheres (with the radius R) arranged in the space with the cubic periodicity. To make a 3D analogy with Model C described in the previous sections, we must find a matrix transform **F**(*x*, *y*, *z*) whose columns, denoted u, v, and w, are complex vectorial functions of space, mutually orthonormal and continuous at all points, where u is the normal vector at each point of the sphere's surface and v and w are tangential.

As the first step we choose the distribution of these vectors on the sphere. Defining the spherical coordinates

$$\begin{aligned} \begin{array}{c} \begin{array}{c} \begin{array}{c} \begin{array}{c} \end{array} \end{array} \end{array} \end{array} \quad \begin{array}{c} \begin{array}{c} \begin{array}{c} \end{array} \end{array} \end{array} \quad \begin{array}{c} \begin{array}{c} \begin{array}{c} \end{array} \end{array} \end{array} \quad \begin{array}{c} \begin{array}{c} \begin{array}{c} \end{array} \end{array} \end{} \begin{array}{c} \begin{array}{c} \begin{array}{c} \end{array} \end{array} \end{} \end{} \end{,} \end{\begin} \end{\begin{array}{c} \begin{array}{c} \begin{array}{c} \end{array} \end{array} \end{aligned} \end{aligned} \end{(86)}$$

$$y = r \sin \theta \sin \phi,\tag{87}$$

$$z = r \cos \theta \tag{88}$$

and the corresponding unit vectors ˆr, ϑ ˆ , and φ ˆ (pointing along the increase of the corresponding coordinate) helps us to define the vectors

$$u\_R(\theta,\phi) = u(R,\theta,\phi) = e^{i(\theta+\phi)}\hat{r},\tag{89}$$

$$v\_R(\vartheta, \phi) = v(R, \vartheta, \phi) = e^{i(\vartheta + \phi)} \frac{1}{\sqrt{1 + \cos^2 \theta}} (\vartheta + i\hat{\phi} \cos \theta), \tag{90}$$

$$w\_R(\theta,\phi) = w(R,\theta,\phi) = e^{i(\theta+\phi)} \frac{1}{\sqrt{1+\cos^2\theta}} (\hat{\theta}\cos\theta - i\hat{\phi}).\tag{91}$$

Obviously, u is normal on the whole sphere's surface. The vector v corresponds to the left circular polarization at the sphere's poles (*z* = ±*R*), to the vertical linear polarization on the sphere's equator (*z* = 0), and is elliptical and continuous in the intermediate ranges. The vector w is simply chosen orthogonal to both u and v. Both v and w are tangential to the sphere's surface.

If we extend the radial dependence of the vectors defined by Equations 89–91 to the entire cubic cell (formally replacing *R* by *r*), then we obtain some analogy with Model B. The obtained matrix **F** has no concurrent discontinuities with *ε*, but it is discontinuous at the center and on the boundaries of the cubic cell. To remove the discontinuity at the cell's center, we will proceed differently.

As the second step, we define the basis vectors at the center of the sphere,

$$u\_0 = u(0, \vartheta, \phi) = \frac{1}{\sqrt{2}}(\hat{x} - i\hat{y}),\tag{92}$$

$$v\_0 = v(0, \theta, \phi) = \hat{\mathbf{z}},\tag{93}$$

$$w\_0 = w(0, \theta, \phi) = \frac{1}{\sqrt{2}}(\hat{x} + i\hat{y}),\tag{94}$$

which are again mutually orthogonal. Then we extend u onto the whole cubic cell,

$$u\_{\rm ext}(r,\theta,\phi) = \begin{cases} \frac{1}{2}(u\_0 + u\_R) + \frac{1}{2}(u\_0 - u\_R)\cos\frac{\pi r}{R} & (r \le R) \\\frac{1}{2}(u\_0 + u\_R) + \frac{1}{2}(u\_0 - u\_R)\cos\frac{\pi(r+D-2R)}{D-R} & (R < r \le D) \end{cases} \tag{95}$$

where *D*(*ϑ*, *φ*) is the distance from the cell's center to its boundary along the ray determined by the spherical angles. The desired unit vector of the polarization basis is then

$$u = A\_u u\_{\text{ext}\nu} \tag{96}$$

where *A*u =[ 1 2 (1 + cos 2 *πr R* )] −1/2 for *r* < *R* and *A*<sup>u</sup> =[ 1 2 (1 + cos 2 *π*(*r*+*D*−2*R*) *D*−*R* )] −1/2 for *r* > *R* is a scalar function ensuring that u becomes a unit vector everywhere. We could extend v and w in a similar way,

$$v\_{\text{ext}}(r,\theta,\phi) = \begin{cases} \frac{1}{2}(\boldsymbol{v}\_{0} + \boldsymbol{v}\_{R}) + \frac{1}{2}(\boldsymbol{v}\_{0} - \boldsymbol{v}\_{R})\cos\frac{\pi r}{R} & (r \le R) \\\frac{1}{2}(\boldsymbol{v}\_{0} + \boldsymbol{v}\_{R}) + \frac{1}{2}(\boldsymbol{v}\_{0} - \boldsymbol{v}\_{R})\cos\frac{\pi(r + D - 2R)}{D - R} & (R < r \le D) \end{cases} \tag{97}$$

$$w\_{\text{ext}}(r,\theta,\phi) = \begin{cases} \frac{1}{2}(w\_0 + w\_R) + \frac{1}{2}(w\_0 - w\_R)\cos\frac{\pi r}{R} & (r \le R) \\ \frac{1}{2}(w\_0 + w\_R) + \frac{1}{2}(w\_0 - w\_R)\cos\frac{\pi(r + D - 2R)}{D - R} & (R < r \le D) \end{cases} \tag{98}$$

but it is not clear whether vext or wext are both perpendicular to u and mutually. To ensure it we define

$$v = A\_v(\mathbf{1} - \mathbf{P}\_u)v\_{\text{ext}\prime} \tag{99}$$

where **P**u = uu † is the projector into the space of vectors proportional to u, so that **1** − **P**<sup>u</sup> is the projector to its orthogonal complement. Similarly,

$$w = A\_w(\mathbf{1} - \mathbf{P}\_u - \mathbf{P}\_v)w\_{\text{ext}\prime} \tag{100}$$

where **1** − **P**<sup>u</sup> − **P**<sup>v</sup> is the projector into the space of vectors perpendicular to both u and v. Here *A*v and *A*w are analogously chosen scalar functions ensuring the unit size of the corresponding vectors.

Finally, the matrix of permittivity in the reciprocal space with correct Fourier factorization, directly applicable to the wave Equation 1, becomes

$$\begin{array}{c} \begin{bmatrix} \mathbf{[e]} \end{bmatrix}\_{\mathbf{C}} = \begin{bmatrix} \mathbf{[F]} \end{bmatrix} \begin{bmatrix} \begin{bmatrix} \frac{1}{\varepsilon} \end{bmatrix} \end{bmatrix}^{-1} \begin{bmatrix} \begin{bmatrix} 0 \end{bmatrix} \begin{bmatrix} \mathbf{[0]} \end{bmatrix} \\\begin{bmatrix} \begin{bmatrix} \varepsilon \end{bmatrix} \end{bmatrix} \begin{bmatrix} \mathbf{[0]} \end{bmatrix} \end{array} \begin{bmatrix} \begin{bmatrix} \mathbf{0} \end{bmatrix} \end{bmatrix} \\\begin{bmatrix} \begin{bmatrix} \mathbf{0} \end{bmatrix} \end{bmatrix} \begin{bmatrix} \mathbf{[0]} \end{bmatrix} \begin{bmatrix} \mathbf{[0]} \end{bmatrix} \end{array} \begin{array}{c} \begin{bmatrix} \mathbf{F}^{-1} \end{bmatrix} \end{array} \tag{101}$$

because the first element on the diagonal corresponds to the normal (*u*) components of the fields and the other two correspond to the tangential components (*v*, *w*) of the fields.

# **8. Conclusion**

We have derived the methodology how to apply the Fourier factorization rules of Li (1996) to various photonic crystals. For the case of 1D crystals there is clear consistency of the correct rules with the classical theory of Yeh et al. (1977). For 2D crystals the convergence properties strongly depend on the chosen distribution of the polarization basis; we have shown that it is desirable to choose a distribution as smooth as possible. The method is also usable for periodic elements of any shape, where complicated shapes require complicated distributions of polarization bases. We can also use it to simulate 2D periodic elements made of anisotropic materials, as well as 3D periodic crystals. Moreover, the method can also be used to photonic devices such as photonic crystal waveguides by applying the demonstrated methodology to the device supercell. It is particularly advantageous for devices where high accuracy is required, e.g., for analyzing defect modes near photonic band edges [Dossou et al. (2009); Mahmoodian et al. (2009)], and for large devices for which the available computer memory enables calculations with only a few Fourier components (photonic crystals fibers with large cladding, or asymmetric 3D crystals).

### **9. Acknowledgments**

This work is part of the research plan MSM 0021620834 financed by the Ministry of Education of the Czech Republic and was supported by a Marie Curie International Reintegration Grant (no. 224944) within the 7th European Community Framework Programme and by the Grant Agency of the Czech Republic (no. 202/09/P355 and P204/10/P346).

### **10. References**



**Photonic Crystals - Introduction, Applications and Theory** Edited by Dr. Alessandro Massaro

ISBN 978-953-51-0431-5 Hard cover, 344 pages **Publisher** InTech **Published online** 30, March, 2012 **Published in print edition** March, 2012

The first volume of the book concerns the introduction of photonic crystals and applications including design and modeling aspects. Photonic crystals are attractive optical materials for controlling and manipulating the flow of light. In particular, photonic crystals are of great interest for both fundamental and applied research, and the two dimensional ones are beginning to find commercial applications such as optical logic devices, micro electro-mechanical systems (MEMS), sensors. The first commercial products involving twodimensionally periodic photonic crystals are already available in the form of photonic-crystal fibers, which use a microscale structure to confine light with radically different characteristics compared to conventional optical fiber for applications in nonlinear devices and guiding wavelengths. The goal of the first volume is to provide an overview about the listed issues.

# **How to reference**

In order to correctly reference this scholarly work, feel free to copy and paste the following:

Roman Antos and Martin Veis (2012). Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals, Photonic Crystals - Introduction, Applications and Theory, Dr. Alessandro Massaro (Ed.), ISBN: 978-953-51-0431-5, InTech, Available from: http://www.intechopen.com/books/photonic-crystalsintroduction-applications-and-theory/fourier-factorization-in-the-plane-wave-expansion-method-in-modelingphotonic-crystals

# **InTech Europe**

University Campus STeP Ri Slavka Krautzeka 83/A 51000 Rijeka, Croatia Phone: +385 (51) 770 447 Fax: +385 (51) 686 166 www.intechopen.com

# **InTech China**

Unit 405, Office Block, Hotel Equatorial Shanghai No.65, Yan An Road (West), Shanghai, 200040, China Phone: +86-21-62489820 Fax: +86-21-62489821

© 2012 The Author(s). Licensee IntechOpen. This is an open access article distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.